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Abstract 

Stochastic reaction-diffusion models have become an important tool in studying how both noise in 
the chemical reaction process and the spatial movement of molecules influences the behavior of biological 
systems. There are two primary spatially-continuous models that have been used in recent studies; the 
diffusion limited reaction model of Smoluchowski, and a second approach popularized by Doi. Both 
models treat molecules as points undergoing Brownian motion. The former represents chemical reactions 
between two reactants through the use of reactive boundary conditions, with two molecules reacting 
instantly upon reaching a fixed separation (called the reaction-radius). The Doi model uses reaction 
potentials, whereby two molecules react with a fixed probability per unit time, A, when separated by less 
than the reaction radius. In this work we study the rigorous relationship between the two models. For 
the special case of a protein diffusing to a fixed DNA binding site, we prove that the solution to the Doi 
model converges to the solution of the Smoluchowski model as A — > oo, with a rigorous 0(A~ 5 +e ) error 
bound (for any fixed e > 0). We investigate by numerical simulation, for biologically relevant parameter 
values, the difference between the solutions and associated reaction time statistics of the two models. As 
the reaction-radius is decreased, for sufficiently large but fixed values of A, these differences are found to 
increase like the inverse of the binding radius. 

1 Introduction 

Stochastic reaction-diffusion models have become a popular tool for modeling biological systems in which 
both noise in the chemical reaction process and the spatial diffusion of molecules play an important roll. Such 
models have been used in a multitude of recent studies, examining the dynamics of synaptic transmission 33 ; 
the MinCDE system in bacteria [TO] ; how proteins search for DNA binding sites [33] ; and studies of signaling 
in the cell membrane [5]. There are three primary stochastic reaction-diffusion models that these studies 
have made use of; the diffusion limited reaction model of Smoluchowski [351 HI]; what we call the Doi 
model [38l[5j|6], and the reaction diffusion master equation (RDME) [Tol \12\ [27 ] . 

In the Doi and Smoluchowski models, the positions of molecules are represented as points undergoing Brown- 
ian motion. Bimolecular reactions between two molecules in the Doi model occur with a fixed probability per 
unit time when two reactants are separated by less than some specified "reaction radius" . The Smoluchowski 
model differs by representing bimolecular reactions in one of two ways; either occurring instantaneously, or 
with fixed probability per unit time, when two reactants' separation is exactly the reaction- radius [35(I28|. In 
this work we focus on the former case (often called a pure absorption reaction). For both models unimolcc- 
ular reactions represent internal processes. They are assumed to occur with exponentially distributed times 

* Department of Mathematics and Statistics, Boston University, 111 Cummington St., Boston, MA 02215 (ag- 
banusi@math.bu.edu) 

^Department of Mathematics and Statistics, Boston University, 111 Cummington St., Boston, MA 02215 (isaac- 
son@math.bu.edu) 



1 



2 



based on a specified reaction-rate constant. For general chemical systems, both the Doi and Smoluchowski 
models can be described by, possibly infinite, systems of partial integral differential equations (PIDEs) for 
the probability densities of having a given number of each chemical species and the corresponding locations 
of each molecule. 

The RDME is spatially discrete, and given by a, possibly infinite, system of ODEs for the numbers of each 
chemical species located at each lattice site. It can be interpreted as an extension of the non-spatial chemical 
master equation (CME) [TH] G21 H51 12S] model for stochastic chemical kinetics. Formally, the RDME has been 
shown to be an approximation to both the Doi model [22] and the Smoluchowski model [221 HS1 H31 120] for 
appropriately chosen lattice spacings. These approximations break down for systems involving bimolecular 
reactions, in two or more dimensions, in the limit that the lattice spacing approaches zero |23LI20| . (Recently 
we have suggested a new convergent RDME (CRDME) that converges to the Doi model as the lattice spacing 
is taken to zero |24j.) 

There are a plethora of numerical methods and simulation packages that have been developed to study biolog- 
ical systems based on one of the Smoluchowski, Doi, or RDME models. These include the X-p method [TTll31| . 
the CRDME [24], the FPKMC [3[36], MCELL [29], MesoRD [14], Smoldyn fl], STEPS [21], and URDME g]. 
While these numerical methods and software packages have been used in many modeling efforts, it is still 
an open question how exactly the underlying mathematical models they approximate are rigorously related. 
Moreover, to compare the results of modeling studies it would be helpful to understand how to choose pa- 
rameters in the Doi (resp. Smoluchowski) model to accurately approximate the Smoluchowski (resp. Doi) 
model. 

To address these questions, we investigate the rigorous relationship between the Doi and (pure absorption) 
Smoluchowski models. We begin in the next section by giving a brief overview of the general Doi and 
Smoluchowski models for the bimolecular reaction A + B — » C. (The two models are identical for systems 
with only first or zcroth order reactions.) We then consider the special case of a single protein searching for 
a DNA binding site within the nucleus of a eukaryotic cell in Section [3] The binding site is located at the 
origin, and the nucleus is modeled as a concentric sphere of radius R. With the further assumption that 
the initial distribution of the protein's position is rotationally invariant, the Doi and Smoluchowski models 
can be simplified to spherically-symmetric diffusion equations. For the Doi model the diffusing protein is 
allowed to bind with probability per unit time A when within a reaction-radius, rt,, of the binding site. This 
leads to a "reaction potential" in the PDE for the Doi model. In the Smoluchowski model the protein reacts 
instantly upon reaching a separation r\> from the binding site. This leads to the replacement of the reaction 
potential with a zero Dirichlet boundary condition on the sphere of radius around the origin. 

Denote by p\ (r, t) the spherically symmetric probability density that the diffusing molecule is r from the 
origin at time t in the Doi model, and by p(r, t) the corresponding probability density in the Smoluchowski 
model. In Section [4] we numerically calculate the difference between p\(r,t) and p{r,t) for biologically 
relevant parameter values. We find that p\(r, t) — > p(r, t) uniformly in r and t as A —¥ oo, with an empirical 
convergence rate that is 0(A _ 5). The same empirical convergence rate is observed for the corresponding 
binding time distributions and mean binding times of the two models. For sufficiently large, but fixed, values 
of A, as rj, is decreased the difference between the probability densities, binding time distributions, and mean 
binding times increase like rZ . 

These results motivate our studies in Section [5j where we rigorously prove the convergence of p\(r,t) to 
p(r,t) as A — > oo with an 0(A _ 5+ e ) error bound (for all e > 0). Let p n denote the nth eigenvalue of the 



generator of the Doi model (see (15)), with a n the nth eigenvalue of the generator of the Smoluchowski 
model (see (11)). Our approach is to first prove that for A sufficiently large, if a n < M(X) with M(X) a 
specified increasing function of A, then 

C 

\Otn ~ fin] < -T— 

A2 

for any fixed e > 0. The precise statement of this result is given in Theorem |5.1[ This theorem is then 
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used to show the corresponding eigcnfunction convergence result in Lemma 5.1 Denote by - 0n( r ) the nth 
eigenfunction of the generator of the Doi model (15), and 4> n (r) the corresponding eigenfunction of the 



Smoluchowski model (11). We prove for A sufficiently large and fj, n < M(A) that 



sup \(j) n (r) 

re[r b ,R] 



Mr)\ < 



c 



Finally, the preceding two results are used to prove that 



sup sup \p\(r,t) - p{r,t)\ < -j 

t£(5,oo) re(r b ,-R) A^ 



C 



(1) 



for any e > and 5 > fixed. The precise statement of this result is given in Theorem 5.2 



It should be noted that the approximation of Dirichlet boundary conditions by reaction potentials is a well- 
studied problem in the context of the large-coupling limit in quantum physics [19j . Our approach of proving 
convergence through successive eigenvalue, eigenfunction, and PDE solution estimates differs from the more 
standard resolvent and path integral estimates 37, 3, 2J. A similar convergence rate of the Doi model solution 
to the Smoluchowski model solution was proven in L 2 (M. 3 ) in [1] (as opposed to the uniform convergence 
rate we prove in a spherical domain). Denote by l[o jrb ]( r ) the indicator function of the interval, [0, rv,]. 
Convergence rates for eigenvalues of the general one-dimensional operator —j^i + V(x) + XW(x) as A — > oo, 



for x £ K, were proven in 
arising in the Doi model, —-^pi 
directly. 



|17j . In contrast, in Section [5] we study the spherically symmetric operator 
Al[o jrb ](V) for r £ [0, R) with a Neumann boundary condition at R, 



2 d_ 

r dr 



2 General Doi and Smoluchowski Models 

We first illustrate how the bimolecular reaction A + B — > C would be described by the Doi and Smoluchowski 
models in R 3 . In the Doi model, bimolecular reactions are characterized by two parameters; the separation 
at which molecules may begin to react, rb, and the probability per unit time the molecules react when within 
this separation, A. In the Doi and Smoluchowski models, when a molecule of species A and a molecule of 
species B react we assume the C molecule they produce is placed midway between them. 

We now formulate the Doi model as an infinite coupled system of PIDEs. Let qf <G M 3 denote the position 
of the Ith molecule of species A when the total number of molecules of species A is a. The state vector of 
the species A molecules is then given by q a = (qf, . . . , q") G M. 3a . Define q b and q c similarly. We denote by 
f(a,b,c) (^qa^ qb qc £j foe probability density for there to be a molecules of species A, b molecules of species B, 
and c molecules of species C at time t located at the positions q a , q h , and q c . Molecules of the same species 
are assumed indistinguishable. In the Doi model the evolution of /( a ' fc > c ) is given by 

- ± W - (q a ,q\ Q c , t) = (L + R) (q a , q\ q c , t) . (2) 

Note, with the subsequent definitions of the operators L and R this will give a coupled system of PIDEs 
over all possible values of (a, b, c). More general systems that allow unbounded production of certain species 
would result in an infinite number of coupled PIDEs. The diffusion operator, L, is defined by 

(a b c \ 

d a a? + d b < + d c A n) f {aAc) (q a ,q b , q c , t) , (3) 
1=1 m=l n=l ) 

where A° denotes the Laplacian in the coordinate qf and D A the diffusion constant of species A. D B , D c , 
A^, and A° are defined similarly. 
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To define the reaction operator, R, we must introduce notations to represent removing or adding a specific 
molecule to the state q a . Let 

q a \ qf = (<??, . . . , q?_i, Ql+1, ■ ■ •,<£), q a U q = (</?, ...,q a a ,q). 

q a \q will denote q a with any one component with the value q removed. Denote by l[ r < rb ](r) the indicator 
function of the interval [0, ft,], and let Bf = {q £ M. 3 \ \q — qf\ < label the set of points a reactant could 

be at to produce a molecule of species C at qf. The Doi reaction operator, R, is then 



(Rf^)(q a ,q h ,q c ,t)=\ 



E / fi»+hb + l,c-l) (<f y g gb {J {2g c _ qhq c\ ql t ) dB c 

-EE W] (Itf - ^ |)/ (a ' 6 ' c V, q\ q c , t) 



(4) 



i=i i'=i 

In the Fock space representation of [HE], the term i?/( a > b ' c ) is called a reaction potential [6]. 

The Smoluchowski model would modify ^ by adding a reactive boundary condition and changing the 
reaction operator Q . The reactive boundary condition, often called a "pure absorption" boundary condition, 
is the Dirichlet boundary condition that 

/( a,M ( g « >g 6 )g c jt ) = o, if \qf -q b m \ = r h , for any l,m. (5) 

We dehne by Sf = {q £ M. 3 \ \q — q^\ = the set of points reactants could be at to produce a molecule 

of species C at qf. For q £ Sf, denote by j( a+1 ' b + 1 - c - 1 ) ( g a Uq,q b (j (2qf -q),q c \ qf, t) the diffusive flux 
along the inward normal to the surface Sf from f( a + 1 > b + 1 - c - 1 ) (q a (J q,q b U (2qf — q) ,q c \ qf, t) , immediately 
prior to a reaction producing qf. The reaction-operator, R, is then 

(jR/ (a,M } ^ qb g c t )=f2 [ j(«+l,b+l,c-l) (g a u g gb y {2g c _ qhq c\ q c Q 

z=i -/^e-sp 

Note, in the Smoluchowski model if the molecules are instead assumed to react with some fixed probability 
per unit time upon collision, a "partial absorption" Robin boundary condition would be used, see [281 113) . 
We do not consider this more general reaction mechanism in this work. 



3 Diffusion of a Protein to a Fixed DNA Binding Site 

To study the rigorous relationship between the Doi and Smoluchowski models we investigate the special case 
of the chemical reaction A + B — > 0, with only one molecule of species A and one molecule of species B. We 
further assume the molecule of species B is located at the origin and stationary (D B = 0). The molecule 
of species A is assumed to move within a sphere of radius R centered on the B molecule. We denote the 
diffusion constant of the A molecule by D and the reaction radius by r^. While idealized, this special case 
can be interpreted as a model for the diffusion of a DNA binding protein to a fixed DNA binding site (located 
at the center of a nucleus). 

We now formulate the Doi and Smoluchowski models in this special case, with the additional assumption 
of spherical symmetry. This assumption will hold whenever the initial distribution of the A molecule is 
spherically-symmetric about the origin. Denote by p{r, t) the spherically-symmetric probability density that 
the A molecule is a distance r from the origin and has not reacted with the B molecule at time t. Then the 
Smoluchowski model ([Si), ^ reduces to 

% = DA rP , r h <r<R,t>0, (7) 
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where A r denotes the spherically symmetric Laplacian in three-dimensions, 



r 2 dr 



A r = ~~ [ r 2 



This equation is coupled with the reactive Dirichlet boundary condition, and zero Neumann boundary 
condition (so that the molecule remains in the "nucleus") 



p(r h ,t) = 0, 



dp 
dr 



(R,t)=0, t>0. 



Finally, we denote by g(r) the initial condition, p(r, 0) = g(r), with the normalization J r g(r)r 2 dr = 1. 

Let p\ (r, t) label the corresponding spherically-symmetric probability density for the Doi model. In the 
special case we are considering, the PIDEs for the Doi model ([2|, Q, Q reduce to 



^ = DA r p\ - A l [0iPb] (r)p A (r, t), < r < R, t > 0, 
with the Neumann boundary condition, 



(8) 



dr 



(R,t) = 0, t>0, 



and the initial condition that 



Px(r,0) = g(r) 



g(r), r > r h , 
0, r < r h . 



For simplicity, in what follows, we assume D = 1. Equations ^ and ^ can be solved explicitly by 
separating variables. In solving ^ we impose continuity of the function and its derivative across the surface 
of discontinuity as justified by the results of [18] and [34]. The computations are standard so we give only 
the final results. Denote by (u(r),v(r)) = u{r)v(r)r 2 dr the usual L 2 inner product. We can then write 
the solutions as 



and 



Here 



(r)= 1 
r 



p(r,t) = £a n (s(r), < M'WnWe- Q » t 

n=l 

oo 

71=1 

sin( A /an(i2 - r)) 



(9) 



(10) 



cos(ya^(i? - r)) 



r b < r < i? 



are the eigenfunctions for the Smoluchowski model ([7]), satisfying 

-A r (f> n (r) = a n (f> n (r), r h < r < R, 



(11) 



with the boundary conditions <^ n (n>) = and ^r(-R) = 0. We extend these functions to [0,R] by defining 
<pn(f) = for r € [0, rb]. The corresponding eigenvalues, a n , solve the equation f(a n ) — 0, where 

RJpZ — tsxi(JJi(R — rb)) 



Rp tan( v //i(i? — r b )) + 



(12) 



G 



The eigenfunctions for the Doi model (J8]) are given by 
where 



VCM, < r < r h , 
r h <r<R, 



sin(^(i? - r b )) - cos(^(i? - r b )) / s inh(rVA - p n ) 



r 



sinh(r b% /A - A*n) 
^= sin(^(i? - r b )) - cos(^(i? - r b )) / s i n ( r ^ Mn _ A) 



sin(r b VMn - A) 



/i„ < A, 

Mn > A, 



and 



sin( v //Z^(.R - r)) 

R^/Pn 



cos( v /^(i? - r)) 



They satisfy the equation 

-A r ^> n (r) + Alp^^nW = ^n^n(r), < r < R, 

with the Neumann boundary condition ^r(-R) = 0. 

The Doi eigenvalues, jtx n , solve the equation f(fi n ) — M^n\ A) with 



(13) 

(14) 
(15) 



y A - // tanh (V A ~ /irb ) ' ^ < A ' 
-!— - tan ^ r b) , M > A. 



A(p,X) 



We denote the dependence of the eigenvalues on A by p, n (X). The constants a n and b n are given by 

1 , 1 

n = 



(16) 



4 Difference Between Smoluchowski and Doi Models for Biologi- 
cally Relevant Parameters 

If we interpret ^ and ^ as models for the diffusion of a protein that has just entered the nucleus (ro = R) 
to a DNA binding site, then typically r b would be between .1 and 10 nm [301 El HI , R between 1 and 10 pm, 
and D between 1 and 20^m 2 s _1 . In the following we assume that all spatial units are in micrometers and 
time is in seconds, with R = ro = 1 fim, D = 10/mi 2 s _1 , and A having units of s _1 . We also assume that 
p\(r,0) = p(r, 0) = 6{r — ro)/4irr 2 . p\ and p are then the same as in the previous section, but rescaled by 

(47T)- 1 . 

We numerically evaluated P\(r,t) and p(r,i) in MATLAB using the eigenfunction expansions (10) and d9l). 
The series were truncated at the first term with magnitude smaller than 10~ 10 . In evaluating these series 
numerically it is necessary to calculate a number of the eigenvalues p n and a n . For each term of the eigenfunc- 
tion expansions the transcendental equations for the corresponding Doi eigenvalue, f(p n ) = A(p n ,X), and 
the Smoluchowski eigenvalue, f{a n ) = 0, were solved to 25 digits of precision using the Mathematica Reduce 
function. In Figure la we show the solution to the Smoluchowski model Q, p(r,t), when r b = 10~ 3 pm. 
For short times the solution is localized near R, while a boundary layer develops near r = r b as t increases. 
Figure [Tb| shows the maximum absolute difference between p\ and p for a discrete set of points, 



\\P\{r,t) - p(r,t)||oo = sup sup \p\(ri,tj) - p(ri,tj)\ , 
* j 
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Figure 1: (a) Solution to the Smoluchowski model Q for ?*b = 10 3 /im. Note that the i-axis uses a 
logarithmic scale, (b) Absolute difference in p\(r, t) and p(r,t) as A and are varied. 



where and tj are given in Appendix [5] by Listings [T] and [2] For each fixed value of rb we see that as A — > oo 
the difference between p x and p converges to zero like A^ 1 / 2 . As r h is decreased, the absolute difference 
increases by approximately an order of magnitude for each fixed value of A (for A sufficiently large). 

For many biological models the statistics of the random variable, for the time at which the diffusing molecule 
first binds to the binding site, are of interest. For example, in [26] we studied how this time was influenced 
by volume exclusion due to the spatially varying density of chromatin inside the nucleus of mammalian cells. 
We subsequently denote these random variables by T Doi and T Smol . Statistics of these random variables can 
be calculated from the cumulative distribution functions 



r R 

Prob [T Doi < t] = 1 - 4tt / px(r, t)r 2 dr, 
Jo 



and 



Prob [T Smo i < t] = 1 - An [ p{r, t)r 2 dr. 



We evaluated Prob pboi < t] and Prob [Ts mo i < t] by analytically integrating the eigenfunction expansions ( 10 ) 
and (|9| and evaluating the truncated series in MATLAB. (Using the same method as described above for 
evaluating p\ (r, t) and p(r,t).) 



Figure 2a shows Prob [Tg mo i < t] for varying and demonstrates a constant increase in the binding time as 



rb is decreased (on a logarithmic scale). Figure 2b shows the absolute difference in binding time distribu- 
tions, 

j| Prob [T Doi < t] - Prob [T Smo i < t] \\ x = sup |Prob [T Doi < tj] - Prob [T Smol < , 

tj 

where tj is given by Listing [2] in Appendix |b| We again observe an empirical A~ 5 convergence rate of 
Prob [Tooi < t] to Prob [Tg mc ,i < t] as A —> 00. For a biologically relevant binding radius of 10~ 3 /im, when 
A = 10 11 s _1 the absolute difference between the two distributions is on the order of 10~ 3 . 



In addition to the probability densities and binding time distributions, we also examined the mean time for 
the diffusing protein to find the binding site (when starting a distance rp from the origin) . The mean binding 
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time for the Doi model can be found in a simple closed form by solving the corresponding mean first passage 
time problem (29), see Appendix [A] and is given by 



Epbod (r ) 
Where, for A = X/D, we have that 
u~(r ) 

u+(r ) 



Prob [T Doi < t] dt 



u (r ), r Q < r b , 
u + (r ), r > TV 



(17) 



M 


sinh(\ 


r Xr Q ) 


)\ 


R 3 -rl 


DX \ 


ra 




3D (r h \ 


/Acosh('\/Xrb) 


r l- r o 


R 3 


' 1 


1 


1 

H r - 


| / sinh(v / Ir b ) 


6D 


+ 3D 




''{) 


DX 





i? 3 - r 



3D ^r b VAcosh(VTr b ) - sinh(V%' b ) 
Similarly, the mean binding time in the Smoluchowski model can be found by solving the mean first passage 

(18) 



time problem (30), and is given by 



E Psmol] M 



6D 



i? 3 
3D 



1 1 



Figure 3a shows the mean binding time in the Doi model as A is varied for ro = R. 
The difference between the two mean binding times, for r > r b , is then given by 



|E[T Doi ] (ro) - E [Tsmoi] Ml = 



1 

DX 



sinh(v Ar b ) 



''b 



R 3 



3D ^r b v // Xcosh('\/Xr b ) — sinh(V%"b) 



(19) 



This difference is 0(X a ), consistent with the results we prove in the next section on the uniform convergence 
of p\ to p as A — > oo. (19) also demonstrates that for large, but fixed values of A the absolute difference in 
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mean binding times will increase like r b 1 as r\> is decreased. In Figure 3b the relative difference, 



|E[T S 



SmolJ 



E [T Do 



E[T f 



SmolJ 



(20) 



is graphed as A is varied for r 
10 



R, illustrating the A 1 I 2 convergence of E [Ibm] to E [Tsmoi]- For r^ = 

li „-i 



3 /xm the mean binding times differ by less than 1% when A = 10 Ai s 



Each of Figures lb 2b and 3b illustrate a A 1 / 2 empirical convergence rate, consistent with the bound |T]) 
we prove in the next section (Theorem 5.2 1. Moreover, they demonstrate that as rb is decreased larger values 



of A are required to ensure the absolute difference between the two models remains below a fixed tolerance. 
In all three cases, once A is sufficiently large that the A -1 / 2 convergence rate can be observed, the absolute 
difference appears to scale like r7 l as rb is decreased. 



5 Rigorous Convergence Results 

In this section we study the rigorous relationship between the Doi model Q and Smoluchowski model ([7]). 
In Subsection |5.1| we derive a rigorous error bound on the rate of convergence of the Doi eigenvalues, 



H n {\), to the Smoluchowski eigenvalues, a„, as A — > oo. In Subsection 5.2 we obtain similar estimates 
for the convergence of the Doi eigenfunctions, ipn Ut , to the Smoluchowski eigenf unctions, <j) n . We obtain 
our main result in subsection |5.3| where we use these eigenvalue and eigenfunction estimates to show the 
uniform convergence in space and time of the solution to the Doi model (JsJ , p\(r,t), to the solution of the 
Smoluchowski model ([7]), p(r,t), as A — > oo. (With the error bound ([I]).) 



5.1 Eigenvalue Estimates 

In this subsection we derive estimates for the difference between the Doi, /i n (A), and Smoluchowski, a„, 
eigenvalues. We start by proving some properties of the functions A([i, A) and f(p) that we will find 
useful: 

Proposition 5.1. A(/i, A) is monotone increasing. Furthermore for fi < A, A is positive. 
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Proof. Positivity is trivial since tanh is positive. Note also that .4(0, A) = tanh(\/AYb)/V / A- A simple 
computation shows that for /i < A 

d tanh (^A - /xr b ) - r h y/X - \x (sech 2 (^A - £tr b )) _ sinh (2^/A - ^r b ) - 2r b ^A - /x 

^ (W )= 2(A- M ) 3/2 = 4cosh 2 (yA^r b )(A- M ) 3/2 

and for /i > A 

d 



A( W A) = 



r bV/" — A (sec 2 (\//U - Ar b )) - tan (^M - Ar b ) 2r h y/ fi - A - sin (2-^ - Ar b ) 



4 cos 2 (y/n - Ar b ) (/j - A) 



dM" v "'' v 2(A- M ) 3/2 

The result follows since for u > 0, sinh(u) > w and sin(u) < u. 

Let the vertical asymptotes of /(/i) be denoted by /3 n . They satisfy the equation 

i?/3„tan(y^(i?-r b )) + ^ = 0; /?„ > 0. 
Proposition 5.2. W^e have the following 

1. < ai < /3i < a 2 < . . . < /3„ < a„+i < . . . 

2. /'(/i) < and f(p) > on [0, ai) and (/?„, a n +i) /or n > 1. 

Proof 1. Let k = 1 — r b /i? and a; = We make the change of variable in /(/x) and obtain 

J? (a; - tan(/ta;)) ^ 



3/2 



□ 



N(x) 

x 1 tan(/tx) + x xtan(ra) + 1 



Let d n be such that N(d n ) = and 77„ be such that D(j] n ) = 0. In terms of the old variable, we have 
d n = R^foLn and rj n = R\f(i n ~. N(x) = imply that t&n(nx) — x and D(x) = imply that tan (/ex) = — -. 
Note that the functions tan(/cx), x and — \ are all monotone increasing. Finally if we let 9 n — ^r(2n — 1) 
be the vertical asymptotes of tan(/cx) one easily checks that we have 

< di < 6>i < T}1 < . . . < d n < 6 n < rj n . . . 

This proves 1. 

2. N > on (J (6>„,d„+i) U (0,di), and N < on (J (d n ,6>„). Similarly, D > on |J (t?„,0„ + i) U (OA) 

n>l n>l n>l 

and D < on (6 n ,r] n ). Thus it follows that /(a;) > if and only if x E [J (?/■„, d„ + i) U (0, di). Next, we 

n>l n>l 

show that /' < 0. Note 

i? (tan(«x) — kx sec 2 (ki)) R (sin(2/cx) — 2kx) 



N'(x) = 



and 



D'(x) = tan(«x) + Kxsec 2 (/cx) 



2x 2 cos 2 (kx) 
sin(2/ta;) + 2kx 



2 cos 2 (kx) 

Since |sin(0)| < \9\ it follows that N'(x) < and that D'(x) > 0, with equality in both only when x = 0. 
From what we have shown above, /(a;) > if and only if both N(x) > and D{x) > so that 



/'(/*) = /» ' 



dx 
d/j, 



D{x)N'{x) - N{x)D'{x) 
D(xf 



da; 
d/j, 



< 0. 



□ 
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Proposition 5.3. The Doi eigenvalues /U„(A) satisfy for n such that //„(A) < A 

< /Ui(A) < ai < Pi < fj, 2 (X) < a 2 < ... < (in-i < M«(A) < a n 

Proof This follows from the fact that /(/x) is decreasing wherever it is positive and that A(/i; A) is positive 
and increasing for \i < A. □ 

We will also have need for the following 

Proposition 5.4. Let {7„} denote the eigenvalues for the (positive) radially symmetric Laplacian on [0, R), 
with zero Neumann boundary conditions on the ball of radius R. Let a n and \i n be as above. Then the 
following hold: 

1. The fi n (X) are continuous and monotone increasing in X for all n> 1 

2. For all n>l and any fixed X, we have that 

M0)=7«= (^^) 2 <MA) 

Remark 5.1. The proof is a straightforward application of the variational minimax principle of Poincare. We 
do not show it here. 

Proposition 5.5. For any L e R + , define the index set A(L) = {n\a n < L}. Lf we let \A(L)\ = c&rd(A(L)) 
then given S > there exists constants C^((5) and C^S) such that for L > 5 

Ct{5)y/L < \A(L)\ < C%(5)^L (21) 

Proof. Write L = R\fL and again define k = 1 — r\,/R. Then |-4(X)| is just the number of solutions to 
tan(/ta;) = x which lie in the interval [0, L]. This number is well approximated by the number of vertical 
asymptotes of tan(«a;). It then follows that 

kL . ...... kL 

1 < \A(L)\ < — + 1 

TT TT 

so that 

If L > S the choice C{ (<5) = — — — ^= and (5) = — — — H — ^= gives the proposition. 

7T y/8 71" VS 



□ 



Remark 5.2. In practice we will choose 8 = Att 2 /(R — n,) 2 so that CI — (R — rb)/27r. 



We now give our main convergence estimate for the eigenvalues of the Doi model. The following theorem 
can be regarded as the heart of the subsequent computations. 

Theorem 5.1. Let < cr < \ and define M{\) = K a \ a for Kq > 1. For any fixed a G (0, a ] there exists 
A > such that for A > A 0; JV ^ < \. Then for a n < M(A), fi n (X) —> a n , OiX-^- 2 ^). 

Remark 5.3. Note that in the remainder C will denote an arbitrary constant that may depend on R, r^, and 
Aq. We will also subsequently assume Aq > 1. 
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Proof. Recall that the Doi eigenvalues /i n (A) satisfy 

/(/x„(A))=A( Mn (A),A). 



As before let k — l—r^/R and let x — R^fji. Recalling the definitions of D(x) and N(x) from Proposition 5.2 
define 



B(x; A) = D(x)A(x, A) 



[xtan(«;a;) + 1] tanh ^Vb \J ^ ^ 



It follows that the rescaled Doi eigenvalues x n (X) satisfy 

N(x n (\)) = B(x n (\);\). 

For A'o > 1 we choose A > (^Kq)^ so that < |. Recall f(fi) = f(x). We restrict to {a; : < 

M(A), /(a;) > 0}, and let h(u) = (1 — u)^ 1 / 2 . Since h is monotone it follows that 



Ai? 2 



_ V fl 2 A 

vA 



< 



/A x/A 



As shown in Proposition 5.2 /(x) > implies — A < tan(Kx) < a; so that 



|^;A)1< [1 + Xtan ^ )] ^< [1 + X2] ^ 



Vx 



Write x n (X) — d n — e„. (Recall d n = R^/a^ and r\ n — R^fJJ^.) Then 

N{d n -e n ) = B(x n [X);X). 
Applying the mean value theorem we get, for some e n £ (x n (X), d n ), that 

N(d n ) - N'(e n )e n = B(x n (X);X). 

As N(d n ) = and N'{e n ) < 0, 



|A'(e„)|e n < 



[1 + 4]^ 
\/A 



which gives that 



vj]l + 3%] (2e^cos 2 ( Ke „)) yg [1 + d 2 n ] e n C[l + d 2 n ) 



Ry/X (2ne n - sin(2 Ke „)) " nRy fx (i _ 5^*1) y/\ (l - si °^ e ") ) 



(22) 



For any fixed c € (0, KX\(1)), monotonicity of the eigenvalues implies < < a;i(Ao) < x n (X) < e n < d ni 

and therefore c < 2ne n . Define 1(9) = 1 — ' ' ^ - . It follows easily that for # > c there exists < to < 1 such 
that 

to < Z(0) < 1 



Using this bound in ( 22 1 , in the original unsealed variables we find that 

CR^n [1 + R 2 a n ] 



Vx 



(23) 
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which implies 

2Ca n [1 + R 2 a n ] 

For a n < M(A) = K \ a , M (A) > 1 implies 

C(M(A)) 2 _ CX 2 



0C n ~ < 



□ 



Remark 5.4. Of interest is the possibility of tighter estimates here. In fact, one can show that if f"{fJ.) > 
wherever /(/x) > we actually have 

l^n _ 1 ■ 

A2 



Theorem |5.1| and Proposition |5.4| immediately implies 

Corollary 5.1. For any fixed n, we have that /^ n (A) converges monotonically to a n as A — > oo. 



5.2 Eigenfunction Estimates 

In this subsection we carry over the estimates for the eigenvalues obtained in the last subsection to obtain the 
uniform convergence in r of the eigenfunctions as A — > oo. Though unstated, in the remainder all theorems 



and lemmas include the assumptions of Theorem 5.1 



Lemma 5.1. The (unnormalized) Doi and Smoluchowski eigenfunctions satisfy 



sup |0n(r)-C*(O| 

r£[r b ,R] 



o(a-<*-*>) 



for n n < M(A) = K X a as A -> oo. 



Proof. 



\Mr)-r n ut (r)\ < 



sm(y/a n (R- r)) ^ sin(y^(i? - r)) 



cos(i/jJ^(R — r)) — cos( v / o^"(i? — r)) 



Note that 



7 < — 



:= 7 + 77 



sin(^/a^(i? - r)) sm( v / /^(7i-r)) 



R\fct7i 



R\JOL n 



1 



sin( v //I^(i? - r)) sin( v //^(i?- r)) 



RJohi 



RJJhx 



I a + h- 



We find 



7„ < 



1 



Rr b y/a n 
Similarly we have 



2 sm [R — r) cos (7c — r) — 



< 



RrbJc^n 



h < 



sm(y/JI^(R - r))| 



1 1 



< 



(R - r h ) 
Rr h 
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Combining these with ( 23 1 



For 77 we have 



< C [1 + R 2 a n ] < CM(A) _ CiT 



'A A5 



77 < 



1 



< 



< 



< 



n 

R-r h 



2 sin (fl-r) — v - sin (R-r) — v - — 



r b \/A 
C(M(A)) = _ CK'I 



Va 



\ 1 3ct 
A2 2 



It follows that 



This concludes the proof. 



A2 2 



< 



d(rb,i?,7f ) 



, l_3ff 
A2 2 



(24) 

□ 



We now prove several uniform properties of the eigenfunctions we will use in the next subsection. 
Lemma 5.2. 1. There exist a Aq, C 2 — C2(rb,R), such that for all A > Aq and n £ Z + 



max sup |-0° Mt (r)|, sup |<£„(r)| < C 2 . 
\re[n,R] re[n,R] I 



(25) 



2. Let b n 



j || 2 2 and a„ — \\(f) n \\ 2 ' 2 . Then there exists C3 such that for all A > A 



maxf sup{a„}, sup{6„}) < C 3 . 

^ 71. 71. ' 



(26) 



Proof. 1. We start by defining for z > 0, and r^ < r < R the auxiliary function 



H(z,r) :-- 



1 



sm(sMR - r)) 
R^fz 



3 (v^(7?-r)) 



Note <p n ( r ) = H(a n ,r) and ip^ u (r) = H(pL ni r). Now for z > z > we have that 



\H{z,r)\ 



1 



&m(yfz(R - r)) 



R^ 



cos(y/z(R — r)) 



< 



< 



»"b 
1 



1 



R4~z 
1 



1 



Ry/zo 

By Corollary |5.1| there exists Aq such that for A > A 



=: C 2 (z ). 



Ot-\ 

«i >j"i(A) >Mi(A ) > y > 0. 

Note that we are using the fact that both the Doi and Smoluchowski eigenvalues can be written in non- 
decreasing order. Choosing zq = ot\/2, proves the first part of the lemma. 
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2. To prove the second part start by defining 



h(z) := / (H(z,r)) 2 r 2 dr 



Once again we have that HV^'lll = M^n) an d Il0n||| = h(a n ). A priori we have that h(z) > for all 
z > 0. An explicit computation shows that for z > 0, h(z) is continuous, lim h(z) = (R — r b )/2 > 0, and 

z— foo 

lim h(z) = (i? 3 — r?)/3i? 2 . With the positivity of h{z) on [0, oo), these results imply that A := inf h(z) > 0. 

z— 

ft then follows that ^ 

On = tt^ — m*t = — r < — =: C*3 



||</>„|| 2 /i(a„) : ^ 



and 



< 



1 



1 1 

< 



c 3 



This concludes the proof of the lemma. 



□ 



These results imply that 

Lemma 5.3. There exists C4 such that for n with /i„(A) < a„ < M(A) 

\K ~ a n \ < i_ 4 3CT ■ 
A2 2 

Proof. We start by noting that 



IIC 



tm£ 1 1 2 



Il0n| 



((C^-^dr 



< 



in 0n + ?/>■ 



< 



d(r b ,i?,if ) 



» 1 3(7 

A2 2 



sup 



r dr 



2C\C 2 (R 3 - rg) 
3As~f 



To get the last line, we have used Lemmas |5.1| and |5.2| A direct computation shows that 



^riii-ii^iiii = iicn 



i=sin( v ^'(i?-r b )) -cos(i/]I^(R - r h ))\ rr b 



3inh(r b \A - A*n) 



sinh 2 (ry/ A — /i„) r 2 dr 



< C (sinh(2r b1 /A - //„) - 2r bv /A - /x„) < C 



VA-/i„sinh (r b VA - £t«) 



A* 



(Here we have used that y/X — /i„ > yA — M(A) > \/f ■) Using the preceding bounds we find, that 



I 

n|l2 



1 

n\\2 



lll'Anlll - llV'nllll 
ll^nlllUVnlll 



< 



< 



out II 2 
n II 2 



- 

llV'nlllll^nlll 



I IIC 



o-ui 1 1 2 
2 



II* 



"112 



ll^||l||0n||I 



CCi , 2C 1 C 2 C 2 (i? J -r b 



A3- 



3A2 2 

^ + |c 1 C 2 Ci(ii 3 -rg) 
A 2 o 



The choice C 4 = CC^ + 2CiC 2 C|( J R J - rg)/3 gives the lemma 



□ 



1G 



5.3 An Error Estimate for the Convergence of the Doi to Smoluchowski Model 

We now show the uniform convergence of the Green's function of the radially symmetric Doi PDE pi to the 
Green's function of the radially symmetric Smoluchowski PDE |7| model. The error bound we give shows 
that the convergence of the Doi model to the Smoluchowski model can not be expected to be faster than 
0(A -1 / 2 ) as A oo. 

Theorem 5.2. Fix < do < \ and let the initial condition g(r) — S(r — r )/r 2 , where r G (r^,i?) . For 
t > 5 > 0, there exists a function u(t) and \q such that for all A > Ao we have: 

sup \p(r,t)-p x (r,t)\<^^. (27) 

r£[r b ,R] A 2 

Moreover, for t £ [S, oo), u(t) is uniformly bounded. 



Proof. The main idea is to use the series representation of the solutions to both models to estimate the error. 
There will be a proliferation of constants which we shall repeatedly and unceremoniously denote by C . For 
re (fb, R), we begin by writing: 

Px (r, t) ~ p(r, t)= Y, bnMroWnH^e-^ 1 - a I ^„(r ) ( /» n (r)e^ Q " t 

{n\a n <M{\)} 

+ Y ^V'n(ro)^r t We-' i " (A)t -anV'n(ro)0n(r) e - Q " t 

{n\a n >M(X)} 

:=I + II. 

We deal with the finitely indexed sum, /, first. Define the index set A\ = {n | a n < M(X)}. From now on, 
for simplicity of presentation, we write ijjn for V 1 ""'. Let I = I a + lb + I c + Id, where 

la = E 6 » MM - M'o)) Mr)e-^ W \ h = Y Mn(r ) (Mr) ~ i(r)) e~^ (A)t , 



I c=Y ( bn ~ °") ( / , n( r o)0n(? , )< 



-Mn(A)i 



Id = ^a„0 n (r o )0„(r) ( 



g-/*r.(A)i _ e -a„i 



Recalling that j n denotes the nth eigenvalue of the radically symmetric Laplacian on [0, R) with a zero 
Neumann boundary condition at R (see Proposition 5.4), we find 

„ oo 

|4| < £ |M \Mr)\ |Vn(ro) ~ <Pn(ro)\e-^ < ]T e~^ . 

A A 2 2 -, 

Ax n=l 



Here we have applied Proposition 5.4 Lemma 5.1 and Lemma 5.2 (in particular (24l, (251, and (261). The 
same argument shows 



\h\ < E K\ \Mr) ~ Mr)\ \Mro)\ e""" 



< 



C 



. l_3o 

A 2 2 



E< 



-Int 



and using Lemma [B~3| too we find 



< 



C 

A2 2 



,-7«< 



Finally, we have that 



\h\<Y\ a -\\Mr)\\K{ro)\ 



I _ g-(«n-^n)< 



,-fJ.nt 
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5.1 



For s > 0, 1 — e s < |s| so that, using the same lemmas as before and Theorem 

\Id\ < C? \a n -fx n \t< -j-^ ]T t e~^K 

We now bound the tail of the series, II. First define the index set B\ — {n \ a n > M(X) }. W e now specify 
the choice Kq > 4ir 2 /(R — r^) 2 which guarantees that K a \ a > An 2 /(R — r^) 2 (see Remark 5.2 after the proof 
of Proposition 5.5). Using the uniform bounds on ip n , <f> n , a n , and b n and Proposition |5 . 5| we find 



\n\ < c x Yl ' 

We thus obtain the error estimate that 



c E 



n>Cl VK^\2 



■ ( ( 

i>ii 



-7«* 



\p(r,t)-p x (r,t)\<C 



A2 2 ' 



-7nt 



1 oo 
AS" 2 - Z i 



(28) 



We estimate the terms in ( 28 ) one at a time. First 

nH-K 2 ' 



E e 7,it = E ex p 



i? 2 



< 1 



exp 



x 2 tn 2 



dx = 1 



i? 



while 



n=2 

oo 



n=l 



< 



n=l 

R 2 — 



eMi(Ao) e7r 2 



»,2 



r^ivr 2 
"ft 2 " 



< C. 



Finally, we bound the third term in ( 28 ) : 



E e " 7 "' = E 



i>A2 



< 



exp 

>A*-1 
oo 

exp 

At -2 
1? 



n 2 tir 2 
ar^Tr 2 " 1 



< , erfc 



R 2 

(A-/ 2 - 2)^tt 



7? 



< 



CR 
v/i/rt 



-c\"t 



Combining the preceding estimates we have 



18 



Here we absorbed the maximum of the many constants into C and 



Let 



C(v) = 



t(t) := C 



1 

Ce 
1 



2a 



1 



- 2 



C(a) 



1 



.y/t V ' t^-\ 

For i > (5 > we see that u(t) is uniformly bounded, concluding the proof. 

Remark 5.5. Note that even for t > S > 0, w(t) — > oo as <r — ► because u(t) blows up like (^=^) ■ 

6 Conclusion 

We have shown for the special case of two molecules that may undergo the annihilation reaction, A + B 



□ 



with one of the two molecules stationary, the solution to the Doi model ( 29 1 converges to the solution of the 
Smoluchowski model (30) as A — > oo. A rigorous asymptotic convergence rate that is <3(A _ 2+ e ) ; for all fixed 
e > 0, was proven for the maximum difference between the two models over all r G (r D , R) and t € (<5, oo) (for 
any fixed S > 0). Numerical evaluation of the exact eigenfunction expansions, binding time distributions, 
and mean binding times illustrated this convergence rate, and demonstrated that for sufficiently large fixed 
values of A the difference between the two models scaled like r^ 1 as was decreased. For biologically 
relevant values of r D , such as the reaction-radius for a protein diffusing to a fixed DNA binding site, it was 
found that A should be chosen at least as large as 10 11 s _1 for the mean binding time in the two models to 
differ by less than 1%. 

There are a number of extensions of the current work that would aid in clarifying the rigorous relationship 
between the Smoluchowski and Doi models. Foremost would be the study of more detailed, biologically 
realistic, models in which multiple diffusing and reacting chemical species are present. Such models require 
the introduction of unbinding reactions, C — > A + B, which in the Smoluchowski model require the use 
of unbinding radii (a separation, greater than r D , at which to place the newly created molecules to avoid 
their immediate rebinding). As numerical methods to solve the Doi and Smoluchowski models have been 
used to study biological systems, understanding how parameters in the Doi model should be chosen so as 
to accurately approximate the Smoluchowski model would greatly aid in comparing predictions by these 
different approaches. 
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A Mean Binding Time 

Let E pboi] denote the mean time at which the two molecules in the Doi model ^ first react when initially 
separated by r . E pboi] can be shown to satisfy [13 [32] 



A ro E pboi] (r ) - A l {r<n>} (r)E pboi] M 



D' 



0<r o <R, 



(29) 



with the boundary condition, 

dE [T Do 



dr Q 



l -(R) = 0. 



(Here A = X/D.) The solution to (29) is given by (17) 



The mean time, E [Tg mo i] (r ), at which the two molecules in the Smoluchowski model ([7| first react when 
initially separated by tq can be shown to satisfy [T51 

A ro E [T Smo i] (r ) = ~, r h <r <R, (30) 

with the boundary conditions, 

E[T Smol ](r b ) = 0, 9E f Sm ° l] (fl)=0. 

or 



The solution to (30) is given by (18). 



B Discrete Space and Time Points 

The spatial evaluation points, r"j, are generated in MATLAB by 



rl = rb + rb*[5e— 6 le 
r = [rl; (.01:. 01:1) 


-5 5e-5 le-4 5e-4 le-3 5e- 
' * (R-rb) + rb ] ; 


-3]'; 




Listing 1: r.; points 




The time evaluation points, tj, are j 


generated in MATLAB by 




t = [le-5 le-4 le-3 


.01:. 01:100 101:1:200 200 


10:1000 1200:200:10000]'; 



Listing 2: tj points 
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